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We present numerical results of three-dimensional simulations for the merger of binary neutron stars 
in full general relativity. Hybrid equations of state are adopted to mimic realistic nuclear equations 
of state. In this approach, we divide the equations of state into two parts as P — Pcoid + Pth- fcoid 
is the cold part for which we assign a fitting formula for realistic equations of state of cold nuclear 
matter slightly modifying the formula developed by Haensel and Potekhin. We adopt the SLy and 
FPS equations of state for which the maximum allowed ADM mass of cold and spherical neutron 
stars is « 2.O4M0 and I.SOMq, respectively. Pth denotes the thermal part which is written as 
Pth = (Pth — l)p(£ — Scold), where p, e, Scoid, and Fth are the baryon rest-mass density, total specific 
internal energy, specific internal energy of the cold part, and the adiabatic constant, respectively. 
Simulations are performed for binary neutron stars of the total ADM mass in the range between 
2.4M0 and 2.8Mq with the rest-mass ratio Qm to be in the range 0.9 ^ Qm < 1. It is found that if 
the total ADM mass of the system is larger than a threshold Mthr, a black hole is promptly formed in 
the merger irrespective of the mass ratios. In the other case, the outcome is a hypermassive neutron 
star of a large ellipticity, which results from the large adiabatic index of the realistic equations of 
state adopted. The value of Mthr depends on the equation of state: Mthr ~ 2.7Mq and ~ 2.5Mq 
for the SLy and FPS equations of state, respectively. Gravitational waves are computed in terms 
of a gauge-invariant wave extraction technique. In the formation of the hypermassive neutron star, 
quasiperiodic gravitational waves of a large amplitude and of frequency between 3 and 4 kHz are 
emitted. The estimated emission time scale is 100 ms, after which the hypermassive neutron 
star collapses to a black hole. Because of the long emission time, the effective amplitude may be 
large enough to be detected by advanced laser interferometric gravitational wave detectors if the 
distance to the source is smaller than ~ 100 Mpc. Thermal properties of the outcome formed after 
the merger are also analyzed to approximately estimate the neutrino emission energy. 


PACS numbers: 04.25.Dm, 04.30.-w, 04.40.Dg 

I. INTRODUCTION 

Binary neutron stars Hi inspiral as a result of the 
radiation reaction of gravitational waves, and eventually 
merge. The most optimistic scenario based mainly on a 
recent discovery of binary system PSRJ0737-3039 0 sug¬ 
gests that such mergers may occur approximately once 
per year within a distance of about 50 Mpc 0. Even 
the most conservative scenario predicts an event rate ap¬ 
proximately once per year within a distance of about 100 
Mpc 0. This indicates that the detection rate of gravi¬ 
tational waves by the advanced LIGO will be ~ 40-600 
yr“^. Thus, the merger of binary neutron stars is one of 
the most promising sources for kilometer-size laser inter¬ 
ferometric detectors 0 0 . 

Hydrodynamic simulations employing full general rela¬ 
tivity provide the best approach for studying the merger 
of binary neutron stars. Over the last few years, numer¬ 
ical methods for solving coupled equations of the Ein¬ 
stein and h ydrody namic equations have been developed 
HI^I^ ITnlTlTLIialT^ll3.ll5t and now such simulations 
are feasible with an accuracy high enough for yielding 
scientific results 0 With the current implementation, 
radiation reaction of gravitational waves in the merger 
of binary neutron stars can be taken into account within 


~ 1% error in an appropriate computational setting 0 
This fact illustrates that it is now a robust tool for the de¬ 
tailed theoretical study of astrophysical phenomena and 
gravitational waves emitted. 

So far, all the simulations for the merger of binary neu¬ 
tron stars in full general relativity have been performed 
adopting an ideal equation of state 01000 
(But, see 0 for a simulation in an approximately gen¬ 
eral relativistic gravity.) For making better models of 
the merger which can be used for quantitative compari¬ 
son with observational data, it is necessary to adopt more 
realistic equations of state as the next step. Since the life¬ 
time (from the birth to the merger) of observed binary 
neutron stars is longer than ^ 100 million yrs 0, the 
thermal energy per nucleon in each neutron star will be 
much lower than the Fermi energy of neutrons 00 at 
the onset of the merger. This implies that for modeling 
the binary neutron stars just before the merger, it is ap¬ 
propriate to use cold nuclear equations of state. During 
the merger, shocks will be formed and the kinetic energy 
will be converted to the thermal energy. However, previ¬ 
ous studies have indicated that the shocks are not very 
strong in the merger of binary neutron stars, in contrast 
to those in iron core collapse of massive stars. The rea¬ 
son is that the approaching velocity at the first contact 
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of two neutron stars is much smaller than the orbital ve¬ 
locity and the sound speed of nuclear matter (e.g., fl^l. 
This implies that the pressure and the internal energy 
associated with the finite thermal energy (temperature) 
are not still as large as those of the cold part. From this 
reason, we adopt a hybrid equation of state in which the 
finite-temperature part generated by shocks is added as a 
correction using a simple prescription (see Sec. II B). On 
the other hand, realisti c eq uations of state are assigned 
to the cold part l2ll| . 

The motivation which stimulates us to perform new 
simulations is that the stiffness and adiabatic index of 
the realistic equations of state are quite different from 
those in the T-law equation of state with T = 2 (hereafter 
referred to as the T = 2 equation of state) which has been 
widely adopted so far (e.g., ^ 3 ) can be expected 
that these differences will modify the properties of the 
merger quantitatively as described in the following. 

Since the realistic equations of state are softer than the 
r = 2 one, each neutron star becomes more compact (cf. 
Fig. EJ. This implies that the merger will set in at a 
more compact state which is reached after more energy 
and angular momentum are already dissipated by grav¬ 
itational radiation. Namely, compactness of the system 
at the onset of the merger is larger. This will modify the 
dynamics of the merger, and accordingly, the threshold 
mass for prompt black hole formation (hereafter Mthr) 
will be changed. 

The adiabatic index of the equations of state is also 
different from that for the F = 2 equation of state. This 
will modify the shape of the hypermassive neutron stars 
^, which are formed after the merger in the case that the 
total mass is smaller than Mj^r- Previous Newtonian and 
post Newtonian studies have indicated that 

for smaller adiabatic index of the equations of state, the 
degree of the nonaxial symmetry of the formed neutron 
star becomes smaller. However, if its value is sufficiently 
large, the formed neutron star can be ellipsoidal. As 
a result of this change, the amplitude of gravitational 
waves emitted from the formed neutron star is signifi¬ 
cantly changed. Since the adiabatic index of the realistic 
equations of state is much larger than that of the F = 2 
equation of state for supranuclear density Hi 110,123 , the 
significant modification in the shape of the hypermassive 
neutron stars and in the amplitude of gravitational waves 
emitted from them is expected. 

The paper is organized as follows. In Sec. II A-C, basic 


^ In this paper, we distinguish the stiffness and the magnitude of 
the adiabatic index clearly. We mention “the equation of state 
is softer (stiffer)” when the pressure at a given density is smaller 
(larger) than another. Thus, even if the adiabatic index is larger 
for the supranuclear density, the equation of state may be softer 
in the case that the pressure is smaller. 

^ The hypermassive neutron star is defined as a differentially ro¬ 
tating neutron star for which the total baryon rest-mass is larger 
than the maximum allowed value o f rig idly rotating neutron stars 
for a given equation of state: See l22l for definition. 


equations, gauge conditions, methods for extracting grav¬ 
itational waves, and quantities used in the analysis for nu¬ 
merical results are reviewed. Then, the hybrid equations 
of state adopted in this paper are described in Sec. II D. 
In Sec. Ill, after briefly describing the computational set¬ 
ting and the method for computation of initial condition, 
the numerical results are presented. We pay particular- 
attention to the merger process, the outcome, and gravi¬ 
tational waveforms. Section IV is devoted to a summary. 
Throughout this paper, we adopt the geometrical units 
in which G = c = 1 where G and c are the gravitational 
constant and the speed of light. Latin and Greek in¬ 
dices denote spatial components (ar, y, z) and space-time 
components (t,x,y,z), respectively: r = \/+ y"^ + z'^. 
Sij{= denotes the Kronecker delta. 

II. FORMULATION 
A. Summary of formulation 

Our formulation and numerical scheme for fully general 
relativistic simulations in three spatial dimensions are the 
same as in m , to which the reader may refer for details 
of basic equations and successful numerical results. 

The fundamental variables for the hydrodynamics are 
p\ rest-mass density, e : specific internal energy, P : 
pressure, : four velocity, and 


where subscripts denote x,y and z, and p, 

the spacetime components. The fundamental variables 
for geometry are a: lapse function, /3^: shift vector, 
7 ij: metric in three-dimensional spatial hypersurface, 
7 = = det( 7 y), 7 ^ = e“"^'^ 7 y: conformal three- 

metric, and Kij'. extrinsic curvature. 

For a numerical implementation of the hydrodynamic 
equations, we define a weighted density, a weighted four- 
velocity, and a specific energy defined, respectively, by 


p* = pau^e®'^. 

(2) 

= hui, 

(3) 

P 

e = hau* - 

(4) 


where h = I -|- e -|- Pjp denotes the specific enthalpy. 
General relativistic hydrodynamic equations are written 
into the conservative form for variables p*, and 

p*e, and solved using a high-resolution shock-capturing 
scheme m- In our approach, the transport terms such as 
di{- ■ ■) are computed by an approximate Riemann solver 
with third-order (piecewise parabolic) spatial interpola¬ 
tion with a Roe-type averaging [^. At each time step, 
au* is determined by solving an algebraic equation de¬ 
rived from the normalization = — 1 , and then, the 

primitive variables such as p, e, and u* are updated. An 
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atmosphere of small density p ~ 10® g/cm^ is added uni¬ 
formly outside neutron stars at t = 0, since the vacuum 
is not allowed in the shock-capturing scheme. The inte¬ 
grated mass of the atmosphere is at most 1% of the total 
mass in the present simulation. Furthermore, we add a 
friction term for a matter of low density ^ 10® g/cm^ to 
avoid infall of such atmosphere toward the central region. 
Hence, the effect of the atmosphere for the evolution of 
binary neutron stars is very small. 

The Einstein evolution equations are solved using a 
version of the BSSN formalism following previous papers 
nimiH: We evolve </>, Aij = e 
and the trace of the extrinsic curvature together with 

three auxiliary functions Fi = using an uncon¬ 

strained free evolution code. The latest version of our 
formulation and numerical method is described in |l 2|| . 
The point worthy to note is that the equation for (j) is 
written to a conservative form similar to the continuity 
equation, and solving this improves the accuracy of the 
conservation of the ADM mass and angular momentum 
significantly. 

As the time slicing condition, an approximate maxi¬ 
mal slice (AMS) condition « 0 is adopted following 
previous papers 0. As the spatial gauge condition, we 
adopt a hyperbolic gauge condition as in [Hill. Suc¬ 
cessful numerical results for the merger of binary neutron 
stars in these gauge conditions are presented in In 
the presence of a black hole, the location is determined 
using an apparent horizon finder for which the method is 
described in (H- 

Following previous works, we adopt binary neutron 
stars in quasiequilibrium circular orbits as the initial con¬ 
dition. In computing the quasiequilibrium state, we use 
the so-called conformally flat formalism for the Einstein 
equation [H- A solution in this formalism satisfies the 
constraint equations in general relativity, and hence, it 
can be used for the initial condition. The irrotational 
velocity field is assumed since it is considered to be a 
good approximation for coalescing binary neutron stars 
in nature [M| . The coupled equations of the field and 
hydrostatic equations |H| are solved by a pseudospec- 
tral method developed by Bonazzola, Gourgoulhon, and 
Marck |H|- Detailed numerical calculations have been 
done by Taniguchi and part of the numerical results are 
presented in |34l| . 


B. Extracting gravitational waves 

Gravitational waves are computed in terms of the 
gau ge-invariant Moncrief variables in a flat spacetime 
|3^ as we have been carried out in our series of paper 
(e.g^^l miHl)- The detailed equations are describe 
in |l2ll37| to which the reader may refer. In this method, 
we split the metric in the wave zone into the flat back¬ 
ground and linear perturbation. Then, the linear part 
is decomposed using the tensor spherical harmonics and 
gauge-invariant variables are constructed for each mode 


of eigen values The gauge-invariant variables of 

I > 2 can be regarded as gravitational waves in the wave 
zone, and hence, we focus on such mode. In the merger 
of binary neutron stars of nearly equal mass, the even- 
parity mode of {I, \m\) = (2, 2) is much larger than other 
modes. Thus, in the following, we pay attention only to 
this mode. 

Using the gauge-invariant variables, the luminosity and 
the angular momentum flux of gravitational waves can be 
defined by 


dE 

dt 

dJ 

dt 


r 

3^ 


E[i« 


.E |2 


t^lm I 


\dtR 


O |2 

im \ 


(5) 


l,m 




l,m 


where Rf^ and are the gauge-invariant variables of 
even and odd parities. The total radiated energy and 
angular momentum are obtained by the time integration 
of dE/dt and dJ/dt. 

To search for the characteristic frequencies of gravita¬ 
tional waves, the Fourier spectra are computed by 


Rim{f) = J e^^^RRlm{t)dt, (7) 


where / denotes a frequency of gravitational waves. Us¬ 
ing the Fourier spectrum, the energy power spectrum is 
defined as 

E (/>o)> (8) 

/>2,m>0 

where for m ^ 0, we define 


Rlmif) = \J\Rln.{fW + \Rl-mifW (ru > 0), (9) 

and use \Rim{-f)\ = \Ri7n{f)\ for deriving Eq. 10). 

We also use a quadrupole formula which is described in 
HE HE ■ shown in , a kind of quadrupole for¬ 

mula can provide approximate gravitational waveforms 
from oscillating compact stars. In this paper, the appli¬ 
cability is tested for the merger of binary neutron stars. 

In quadrupole formulas, gravitational waves are com¬ 
puted from 


hij — 


-p k r) I _ p pkl 

^3 2 ^ 


^kl 


dt^ 


( 10 ) 


where fy and = Sij — riirij {rii = x'‘jr) denote a 
tracefree quadrupole moment and a projection tensor. 

In fully general relativistic and dynamical spacetimes, 
there is no unique definition for the quadrupole moment 
lij. Following HE HE I choose the formula as 



( 11 ) 
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Then, using the continuity equation, the first time deriva¬ 
tive is computed as 

iij = j -\-x'^v^)cPx. ( 12 ) 

To compute Iij , we carry out the finite differencing of the 
numerical result for Iij. 

In this paper, we focus only on ^ = 2 mass quadrupole 
modes. Then, the gravitational waveforms are written as 


in the gauge-invariant wave extraction method, and as 
the corresponding variables, 

■^+ = ixx ~ iyyj ( 22 ) 

Ax = (23) 

in the quadrupole formula. These have the unit of length 
and provide the amplitude of a given mode measured by 
an observer located in the most optimistic direction. 


= 


647r 


{.R 22+(1 + cos^ 9) cos( 2 (^) 


+ -R 22-(1 + cos^ 9) sin(2(/?)} 
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i ?20 sin^ 9 


/lx — 


5 

647r - 


647r 

—i? 22 + cos 9 sin(2(/j) 
i?22- COS0COS(2(/?) 


(13) 


(14) 


in the gauge-invariant wave extraction technique, and 
1 


h+ = - 


0)cos(2v5) 


+ ixy{^ + cos^ 9) sin(2(/j) 


I I j Ixx + lyy 
I I ^zz ^ 


sin^ 9 


hy — — 


Ixx ^ lyy cosgsin(2(^) 
-I- Ixy COS 9 cos(2i^) 


(15) 


(16) 


in the quadrupole formula. In Eqs. and m , we use 
the variables defined by 


R22± = 


dE _|_ pE 
-^22 ^ -^2 —2 


V 2 


i ?20 = -^ 20 ^- 


(17) 

(18) 


For the derivation of /i+ and /ix, we assume that the wave 
part of the spatial metric in the wave zone is written as 

dP = dP + r^[(l -I- h+)d9'^ -\- sin^ 0(1 — h+)dtp^ 

-\-2sin9hxd9d(p\, (19) 

and set i?f = 0 and Ixz = lyz = 0 since we assume the 
reflection symmetry with respect to the equatorial plane. 
In the following, we present 


/ 5 

( 20 ) 



( 21 ) 


C. Definitions of quantities and methods for 
calibration 


In numerical simulations, we refer to the total baryon 
rest-mass, the ADM mass, and the angular momentum 
of the system, which are given by 


M, = / p^.d^x, 


(24) 


M =-Tr f 

Jr—yoo 

IGtt \ 


k\2 






J = 


1 

Stt 


(f pA/e^^dS, 

J r —too 


d^x, (25) 


= / 






d^x, (26) 


where dSj = Pdjrd{cos9)dip, ip = —y{dxy + x{dy)y 
ip = e*^, ph = pau*e, Ji = pui, and denotes the 
Ricci scalar with respect to 7 ^. To derive the expressions 
for M and J in the form of volume integral, the Gauss 
law is used. Here, M* is a conserved quantity. We also 
use the notations M*i and M *2 which denote the baryon 
rest-mass of the primary and secondary neutron stars, 
respectively. In terms of them, the baryon rest-mass ratio 
is defined by Qm = M* 2 /M*i(< 1 ). 

In numerical simulation, M and J are computed using 
the volume integral shown in Eqs. Ill and Since 

the computational domain is finite, they are not constant 
and decrease after gravitational waves propagate to the 
outside of the computational domain during time evolu¬ 
tion. Therefore, in the following, they are referred to as 
the ADM mass and the angular momentum computed in 
the finite domain (or simply as M and J, which decrease 
with time). 

The decrease rates of M and J should be equal to the 
emission rates of the energy and the angular momentum 
by gravitational radiation according to the conservation 
law. Denoting the radiated energy and angular momen¬ 
tum from the beginning of the simulation to the time t as 
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AE{t) and AJ(i), the conservation relations are written 
as 


M{t) + AE{t) = Mo, (27) 

J(t) + AJ(t) = Jo, (28) 


where Mg and Jg are the initial values of M and J. We 
check if these conservation laws hold during the simula¬ 
tion. 

Significant violation of the conservation laws indicates 
that the radiation reaction of gravitational waves is not 
taken into account accurately. During the merger of bi¬ 
nary neutron stars, the angular momentum is dissipated 
by several 10%, and thus, the dissipation effect plays an 
important role in the evolution of the system. Therefore, 
it is required to confirm that the radiation reaction is 
computed accurately. 

The violation of the Hamiltonian constraint is locally 
measured by the equation as 


U 


AV' - f+ 27rpHV'^ 

+ \^Rk^\ + |27rpHV’®| 




(29) 


Following we define and monitor a global quantity 
as 


H J pj,pd^x. (30) 

Hereafter, this quantity will be referred to as the averaged 
violation of the Hamiltonian constraint. 


D. Equations of state 

Since the lifetime of binary neutron stars from the birth 
to the merger is longer than ~ 100 million yrs for the ob¬ 
served systems 0, the temperature of each neutron star 
will be very low (^ 10^ K) [13,113 the onset of merger; 
i.e., the thermal energy per nucleon is much smaller than 
the Fermi energy of neutrons. This implies that for mod¬ 
eling the binary neutron stars just before the merger, it 
is appropriate to use cold nuclear equations of state. On 
the other hand, during the merger, shocks will be formed 
and the kinetic energy will be converted to the thermal 
energy to increase the temperature. However, previous 
studies have indicated that the shocks in the merger are 
not strong enough to increase the thermal energy to the 
level as large as the Fermi energy of neutrons, since the 
approaching velocity at the first contact of two neutron 
stars is much smaller than the orbital velocity and the 
sound speed of nuclear matter. This implies that the 


pressure and the internal energy associated with the fi¬ 
nite temperature are not still as large as those of the cold 
part. From this reason, we adopt a hybrid equation of 
state. 

In this equation of state, we write the pressure and the 
specific internal energy in the form 

P = ^cold + ^th, (31) 

^ — ^cold T £th, (^2) 

where Pcoid and £coid are the cold (zero-temperature) 
parts, and are written as functions of p. Pth and £th 
are the thermal (finite-temperature) parts. During the 
simulation, p and e are computed from hydrodynamic 
variables p* and e. Thus, £th is determined by e — Ecoid- 

For the cold parts, we assign realistic equations of 
state for zero-temperature nuclear matter. In this pa¬ 
per, we adopt the SLy [^ and FPS equations of state 
|19l| . These are tabulated as functions of the baryon rest- 
mass density for a wide density range from ~ 10 g/cm^ 
to ^ 10^® g/cm^. To simplify numerical implementation 
for simulation, we make fitting formulae from the tables 
of equations of state, slightly modifying the original ap¬ 
proach proposed in [^. 

In our approach, we first make a fitting formula for 

S^cold 

ecold(p) = [(l+Plff^ +P3P^^)(1 +P5P^®)^^ - 1] 

X f(-PsP + Pw) 

+ Pl2fj'^^ KpbP - Pw)f(-P9P + Pll) 

+ Pl4P^""/(p9P-Pll), (33) 


where 

f(x) = (34) 

-I- 1 

The coefficients pi (i =1-15) denote constants, and are 
listed in Table I. In making the formula, we focus only 
on the density for p > 10^° g/cm^ in this work, since the 
matter of lower density does not play an important role 
in the merger. Then, the pressure is computed from the 
thermodynamic relation in the zero-temperature limit 

dp 

With this approach, the accuracy of the fitting for the 
pressure is not as good as that in j^. However, the 
first law of the thermodynamics is completely satisfied in 
contrast to that in |^. 

In Fig. 1, we compare Pcoid and £coid calculated by the 
fitting formulae (solid curves) with the numerical data 
tabulated (dotted curves) It is found that two results 


® The tables for the SLy and FPS equations of state, which 
were involved in the LORENE library in Meudon group 
I http://www.lorene.obspm.fr I, were implemented by Haensel 
and Zdunik. 
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i 

Pi (SLy) 

Pi (EPS) 

i Pi (SLy) Pi (FPS) 

1 

0.1037 

0.15806 

9 9 X 10® 

9 X 10® 

2 

0.1956 

0.220 

10 4 

5 

3 

39264 

5956.4 

11 0.75 

0.75 

4 

1.9503 

1.633 

12 0.057 

0.0627 

5 

254.83 

170.68 

13 0.138 

0.1387 

6 

1.3823 

1.1056 

14 0.84 

0.56 

7 

-1.234 

-0.703 

15 0.338 

0.308 

8 

1.2 X 10® 

2 X 10"^ 




TABLE I: The values of pi we choose in units of c = G = Mq = 1. 



lO'o 10' 


1012 ;L0'3 lOl'l lO'" 
p(g/cm3) 



lOio 10' 


1012 ;L0'3 lOl'l 10'3 
p(g/cm3) 


(a) 


(b) 


FIG. 1: Pressure and specihc internal energy as functions of baryon rest-mass density p (a) for the SLy and (b) for the EPS 
equations of state. The solid and dotted curves denote the results by htting formulae and numerical data tabulated, respectively. 


agree approximately. The relative error between two is 
within ^ 10% for p > 10^° g/cm^ and ^ 2% for supranu¬ 
clear density with p2 x 10^^ g/cm^. 

In Fig. 121 we show the relations among the ADM mass 
M, the total baryon rest-mass M*, the central density pc, 
and the circumferential radius R for cold and spherical 
neutron stars in the SLy and FPS equations of state. For 
comparison, we present the results for the F = 2 poly- 
tropic equation of state P = Kpp^ which was adopted 
in [^. In the polytropic equations of state, there ex¬ 
ists a degree of freedom for the choice of the polytropic 
constant Kp. Here, for getting approximately the same 
value of the maximum ADM mass for cold and spherical 
neutron stars as that of the realistic equations of state, 
we set 

Kp = 1.6 X 10® (cgs unit). (36) 

In this case, the maximum ADM mass is about 1.72M0. 
We note that for the F = 2 equation of state, the ADM 
mass M, the circumferential radius R, and the density 
can be rescaled by changing the value of Kp using the 


following rule: 

MocKp^^, RocKp^'^, and pociFp®. (37) 

Hence, the mass and the radius are arbitrarily rescaled 
although the compactness M/R is invariant in the rescal¬ 
ing. 

Figure 121 shows that in the realistic equations of state, 
the central density and the circumferential radius are in 
a narrow range for the ADM mass between ~ O.SMq 
and ~ 1.5Mq. Also, it is found that neutron stars in the 
realistic equations of state are more compact than those 
in the F = 2 polytropic equation of state for a given mass. 
Namely, the realistic equations of state are softer than 
the F = 2 one. On the other hand, the adiabatic index 
din P/d In p for the realistic equations of state is much 
larger than 2 for the supranuclear density [Uliaill. 
These properties result in quantitatively different results 
in the merger of two neutron stars from those found in 
the previous work da 

The thermal part of the pressure Pth is related to the 
specific thermal energy £th = s — EcoW as 

Pth — (Tth l)pS^th; 


(38) 
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(a) (b) 

FIG. 2: (a) ADM mass (solid curves) and total baryon rest-mass (dotted curves) as functions of central baryon rest-mass 

density pc and (b) relation between the circumferential radius and the ADM mass for cold and spherical neutron stars in 
equilibrium. ’FPS’ and ’SLy’ denote the sequences for the FPS and SLy equations of state, respectively. The relations for the 
P = 2 polytropic equation of state with Kp = 1.6 x 10® in the cgs unit are also drawn by the dashed curves. 


where Fth is an adiabatic constant. As a default, we 
set Fth = 2 taking into account the fact that the equa¬ 
tions of state for high-density nuclear matter are fairly 
stiff. (We note that for the ideal nonrelativistic Fermi 
gas, Fth ~ 5/3 li^l. For the nuclear matter, it is reason¬ 
able to consider that it is much larger than this value.) 
To investigate the dependence of the numerical results on 
the value, we also choose Fth = 1-3 and 1.65. The ther¬ 
mal part of the pressure plays an important role when 
shocks are formed during the evolution. For the smaller 
value of Fth — 1, local conversion rate of the kinetic energy 
to the thermal energy at the shocks should be smaller. 


III. NUMERICAL RESULTS 

A. Initial condition and computational setting 

Several quantities that characterize irrotational binary 
neutron stars in quasiequilibrium circular orbits used as 
initial conditions for the present simulations are summa¬ 
rized in Table II. We choose binaries of an orbital sepa¬ 
ration which is slightly larger than that for an innermost 
orbit. Here, the innermost orbit is defined as a close or¬ 
bit for which Lagr ang e points appear at the inner edge of 
neutron stars [331 l4l| . If the orbital separation becomes 
smaller than that of the innermost orbit, mass transfer 
sets in and dumbbell-like structure will be formed. Until 
the innermost orbit is reached, the circular orbit is stable, 
and hence, the innermost stable circular orbit does not 
exist outside the innermost orbit for the present cases. 
However, we should note that the innermost stable cir¬ 
cular orbit seems to be very close to the innermost orbit 


since the decrease rates of the energy and the angular 
momentum as functions of the orbital separation are very 
small near the innermost orbit. 

The ADM mass of each neutron star, when it is in 
isolation (i.e., when the orbital separation is infinity), is 
chosen in the range between 1.2Mq and 1.45M0. Mod¬ 
els SLyl212, SLyl313, SLyl35135, SLyl414, FPS1212, 
FPS125125, FPS1313, and FPS1414 are equal-mass bi¬ 
naries, and SLyl25135 and SLyl35145 are unequal-mass 
ones. For the unequal-mass case, the mass ratio Qm is 
chosen to be ^ 0.9 since all the observed binary neutron 
stars for which each mass is determined accurately in¬ 
deed have such mass ratio Q . Mass of each neutron star 
in model SLyl25135 is approximately the same as that of 
PSRJ0737-3039 0|, while the mass in model SLyl35145 
is similar to that of PSRB1913-I-16 0. The total baryon 
rest-mass for models SLyl313 and SLyl25135 and for 
models SLyl414 and SLyl35145 are approximately iden¬ 
tical, respectively. For all these binaries, the orbital pe¬ 
riod of the initial condition is about 2 ms. This im¬ 
plies that the frequency of emitted gravitational waves is 
about 1 kHz. 

The simulations were performed using a fixed uniform 
grid and assuming reflection symmetry with respect to 
the equatorial plane (here, the equatorial plane is cho¬ 
sen to be the orbital plane). The detailed simulations 
were performed with the SLy equation of state. In this 
equation of state, the used grid size is (633, 633, 317) or 
(377, 377, 189) for {x, y, z). In the FPS equation of state, 
simulations were performed with the (377, 377, 189) grid 
size to save the computational time. The grid covers the 
region —L < x < L, —L < y < L, and 0 < z < L where 
L is a constant. The grid spacing is determined from the 





Model 

Each ADM mass 

pmax 

Qm 

M* 

Mo 

go 

Po 

Co 

Q* 

/o 

SLyl212 

1.20, 1.20 

8.03, 8.03 

1.00 

2.605 

2.373 

0.946 

2.218 

0.103 

1.075 

0.902 

SLyl313 

1.30, 1.30 

8.57, 8.57 

1.00 

2.847 

2.568 

0.922 

2.110 

0.112 

1.175 

0.948 

SLyl35135 

1.35, 1.35 

8.86, 8.86 

1.00 

2.969 

2.666 

0.913 

2.083 

0.116 

1.225 

0.960 

SLyl414 

1.40, 1.40 

9.16, 9.16 

1.00 

3.093 

2.763 

0.902 

2.012 

0.122 

1.277 0.994 

SLyl25135 

1.25, 1.35 

8.29, 8.86 

0.9179 

2.847 

2.568 

0.921 

2.110 

0.112 

1.175 

0.948 

SLyl35145 

1.35, 1.45 

8.85, 9.48 

0.9226 

3.094 

2.763 

0.901 

2.013 

0.122 

1.277 0.994 

FPS1212 

1.20, 1.20 

9.93, 9.93 

1.00 

2.624 

2.371 

0.925 

1.980 

0.111 

1.251 

1.010 

FPS125125 

1.25, 1.25 

10.34, 10.34 

1.00 

2.746 

2.469 

0.914 

1.935 

0.116 

1.309 

1.034 

FPS1313 

1.30, 1.30 

10.79, 10.79 

1.00 

2.869 

2.566 

0.903 

1.882 

0.121 

1.368 

1.063 

FPS1414 

1.40, 1.40 

11.76, 11.76 

1.00 

3.120 

2.760 

0.882 

1.750 

0.134 

1.487 

1.143 


TABLE II: A list of several quantities for quasiequilibrium initial data. The ADM mass of each star when they are in isolation, 
the maximum density for each star, the baryon rest-mass ratio Qm = the total baryon rest-mass, the total ADM 

mass Mo, nondimensional spin parameter qo = Jo/Mg , orbital period Pq, the orbital compactness [Co = the ratio of 

the total baryon rest-mass to the maximum allowed mass for a spherical and cold star (Q* = M ,and the frequency 
of gravitational waves /o = 2/Po. The density, mass, period, and wave frequency are shown in units of 10^‘^g/cm®, Mq, ms, 
and kHz, respectively. 


Model 

Fth 

Grid number 

L 

A 

Ao 

/"merger 

Amerger 

Product 

SLyl212b 

2 

(377, 377, 189) 

77.8 

0.414 

333 

3.1 

97 

NS 

SLyl313a 

2 

(633, 633, 317) 

130.8 

0.414 

316 

3.2 

94 

NS 

SLyl313b 

2 

(377, 377, 189) 

77.8 

0.414 

316 

3.2 

94 

NS 

SLyl313c 

1.3 

(377, 377, 189) 

77.8 

0.414 

316 

3.7 

81 

NS BH 

SLyl313d 

1.65 

(377, 377, 189) 

77.8 

0.414 

316 

3.4 

88 

NS 

SLyl35135b 

2 

(377, 377, 189) 

77.8 

0.414 

316 

3.6 

83 

NS ^ BH 

SLyl414a 

2 

(633, 633, 317) 

130.8 

0.414 

302 

— 

— 

BH 

SLyl25135a 

2 

(633, 633, 317) 

130.8 

0.414 

316 

3.2 

94 

NS 

SLyl35145a 

2 

(633, 633, 317) 

130.8 

0.414 

302 

— 

— 

BH 

FPS1212b 

2 

(377, 377, 189) 

69.5 

0.370 

297 

3.5 

86 

NS BH 

FPS125125b 

2 

(377, 377, 189) 

69.5 

0.370 

297 



NS ^ BH 

FPS1313b 

2 

(377, 377, 189) 

69.5 

0.370 

282 

— 

— 

BH 

FPS1414b 

2 

(377, 377, 189) 

69.5 

0.370 

262 

— 

— 

BH 


TABLE III: A list of computational setting. Lth, L, A, Ao, and /merger denote the adiabatic index for the thermal part, 
the location of outer boundaries along each axis, the grid spacing, the wave length of gravitational waves at t = 0, and the 
frequency of gravitational waves from the formed hypermassive neutron stars, respectively. Amerger denotes the wave length of 
gravitational waves from the formed hypermassive neutron stars Amerger = c//merger. The length and the frequency are shown 
in units of km and kHz. In the last column, the outcome is shown. NS implies that a hypermassive neutron star is produced and 
remains stable at t ~ 10 ms. NS^BH implies that a hypermassive neutron star is formed first, but as a result of gravitational 
radiation reaction, it collapses to a black hole in t 10 ms. BH implies that a black hole is promptly formed. 


condition that the major diameter of each star is cov¬ 
ered with about 50 grid points initially. We have shown 
that with this grid spacing, a convergent numerical result 
is obtained (l^ . The circumferential radius of spherical 
neutron stars with the SLy and FPS equations of state 
is about 11.6 and 10.7 km for M = 1.4Mq, respectively 
(see Fig. El. Thus, the grid spacing is ~ 0.4 km. 

Accuracy in the computation of gravitational wave¬ 
forms and the radiation reaction depends on the location 
of the outer boundaries if the wavelength. A, is larger 
than L For L ^ 0.4A, the amplitude and the ra¬ 
diation reaction of gravitational waves are significantly 
overestimated [l^ . Due to the restriction of the com¬ 

putational power, it is difficult to take a huge grid size in 


which L is much larger than A. As a consequence of the 
present restricted computational resources, L has to be 
chosen as ^ 0.4Ao where Aq denotes A of the I = m = 2 
mode at t = 0. Hence, the error associated with the 
small value of L is inevitable, and thus, the amplitude 
and radiation reaction of gravitational waves are overes¬ 
timated in the early phase of the simulation. However, 
the typical wavelength of gravitational waves becomes 
shorter and shorter in the late inspiral phase, and hence, 
the accuracy of the wave extraction is improved with the 
evolution of the system. This point will be confirmed in 
Sec. HI. 

The wavelength of quasiperiodic gravitational waves 
emitted from the formed hypermassive neutron star (de- 
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FIG. 3: Snapshots of the density contour curves for p in the equatorial plane for model SLyl414a. The solid contour curves 
are drawn for p — 2 x 10^“^ x i g/cm^ (i = 2 ~ 10) and for 2 x 10^'* x 10“° ®* g/cm° (i = 1 ^ 7). The dotted curves denote 
2 X g/cm°. The number in the upper left-hand side denotes the elapsed time from the beginning of the simulation in 
units of ms. The initial orbital period in this case is 2.012 ms. Vectors indicate the local velocity held and the scale is 

shown in the upper right-hand corner. The thick circle in the last panel of radius r ^ 2 km denotes the location of the apparent 
horizon. 


noted by Amerger) is much shorter than Aq and as large 
as L (see Table III), so that the waveforms in the merger 
stage are computed accurately (within ~ 10% error) in 
the case of neutron star formation irrespective of the 
grid size. We performed simulations for models SLyl313, 
SLyl25135, SLyl414, and SLyl35145 with the two grid 
sizes, and confirmed that this is indeed the case. We 
demonstrate this fact in Sec. Ill by comparing the results 
for models SLyl313a and SLyl313b. From the numeri¬ 
cal results for four models, we have also confirmed that 
the outcome in the merger does not depend on the grid 
size. Thus, when we are interested in the outcome or in 
gravitational waves emitted by the hypermassive neutron 
stars, simulations may be performed in a small grid size 
such as (377, 377, 189). 

With the (633, 633, 317) grid size, about 240 GBytes 
computational memory is required. For the case of the 
hypermassive neutron star formation, the simulations are 
performed for about 30,000 time steps (until t ~ 10 ms) 
and then stopped to save the computational time. The 
computational time for one model in such a simulation 
is about 180 CPU hours using 32 processors on FACOM 
VPP5000 in the data processing center of National As¬ 
tronomical Observatory of Japan (NAOJ). For the case 
of the black hole formation, the simulations crash soon 
after the formation of apparent horizon because of the 
so-called grid stretching around the black hole formation 


region. In this case, the computational time is about 60 
CPU hours for about 10,000 time steps. 


B. Characteristics of the merger 

1. General feature 

In Figs. EHSl we display the snapshots of the density 
contour curves and the velocity vectors in the equato¬ 
rial plane at selected time steps for models SLyl414a, 
SLyl313a, and SLyl25135a, respectively. Figure El dis¬ 
plays the density contour curves and the velocity vec¬ 
tors in the y = 0 plane at a late time for SLyl414a and 
SLyl313a. Figures 0 and Ela) indicate typical evolution 
of the density contour curves in the case of prompt black 
hole formation. On the other hand, Figs. a a and 
Elb) show those in the formation of hypermassive neu¬ 
tron stars. 

Figure a displays the evolution of the maximum value 
of p (hereafter Pmax) and the central value of a (here¬ 
after Oc) for models SLyl414a, SLyl35145a, SLyl35135b, 
SLyl313a, SLyl25135a, and SLyl212b (Fig. [7Ka)) and 
for models FPS1414b, FPS1313b, FPS125125b, and 
FPS1212b (Fig. ^b)). This shows that in the prompt 
black hole formation, monotonically decreases toward 
zero. On the other hand, Oc and pmax settle down to cer- 



































10 



X (km) 


X (km) 


X (km) 


FIG. 4: The same as Fig. I^but for model SLyl313a. The initial orbital period is 2.110 ms in this case. 


tain values in the hypermassive neutron star formation. 
For model SLyl35135b, a hypermassive neutron star is 
formed first, but after several quasiradial oscillations of 
high amplitude, it collapses to a black hole due to dissi¬ 
pation of the angular momentum by gravitational radi¬ 
ation. The large oscillation amplitude results from the 
fact that the selfgravity is large enough to deeply shrink 
surmounting the centrifugal force. These indicate that 
the total ADM mass of this model (M « 2.7Mq) is only 
slightly smaller than the threshold value for the prompt 
black hole formation. The quasiradial oscillation of the 
large amplitude induces a characteristic feature in grav¬ 
itational waveforms and the Fourier spectrum (cf. Sec. 

nrm . 

In the case of black hole formation (models SLyl414, 
SLyl35145, SLyl35135, FPS1414, FPS1313, FPS125125, 
and FPS1212), the computation crashed soon after the 
formation of apparent horizons since the region around 
the apparent horizon of the formed black hole was 
stretched significantly and the grid resolution became 
too poor to resolve such region. On the other hand, we 
stopped the simulations for other cases to save the com¬ 


putational time, after the evolution of the formed massive 
neutron stars was followed for a sufficiently long time. At 
the termination of these simulations, the averaged viola¬ 
tion of the Hamiltonian constraint H remains of order 
0.01 (cf. Fig. ITHll . We expect that the simulations could 
be continued for a much longer time than 10 ms if we 
could have sufficient computational time. 

In every model, the binary orbit is stable at t = 0 
and the orbital separation gradually decreases due to the 
radiation reaction of gravitational waves for which the 
emission time scale is longer than the orbital period. If 
the orbital separation becomes sufficiently small, each 
star is elongated by tidal effects. As a result, the attrac¬ 
tion force due to the tidal interaction between two stars 
becomes strong enough to make the orbit unstable to 
merger. The merger starts after about one orbit at t ^ 2 
ms irrespective of models. Since the orbital separation at 
t = 0 is very close to that for a marginally stable orbit, 
a small decrease of the angular momentum and energy 
is sufficient to induce the merger in the present simula¬ 
tions. If the total mass of the system is high enough, a 
black hole is directly formed within about 1 ms after the 
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FIG. 5; The same as Fig. I^but for model SLyl25135a. The initial orbital period is 2.110 ms in this case. 




(a) (b) 


FIG. 6: Snapshots of the density contour curves for p and the local velocity field in the y ~ 0 plane (a) at t = 2.991 

ms for model SLyl414a and (b) at t = 8.621 ms for model SLyl313a. The method for drawing the contour curves and the 
velocity vectors is the same as that in Fig. 
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(a) 


(b) 


FIG. 7: Evolution of the maximum values of p and the central value of a (a) for models SLyl414a (long dashed curves), 
SLyl35145a (dotted long-dashed curves which approximately coincide with long dashed curve for Qc), SLyl35135b (dotted 
curves), SLyl313a (solid curves), SLyl25135a (dotted dashed curves), and SLyl212b (dashed curves), and (b) for models 
FPS1414b (long dashed curves), FPS1313b (solid curves), FPS125125 (dotted-dashed curves), and FPS1212b (dashed curves). 
The dotted horizontal lines denote the central values of the lapse and density of the marginally stable and spherical star in 
equilibrium for given cold equations of state. The reason that the maximum density decreases in the final stage for the black 
hole formation case is as follows: We choose p* as a fundamental variable to be evolved and compute p from p*/(au*e®‘^). In 
the final stage, 4> is very large (> 1) and, hence, a small error in ij) results in a large error in p. Note that the maximum value 
of pt increases monotonically by many orders of magnitude. 



FIG. 8: Evolution of the angular velocity II along x (solid curves) and y (dashed curves) axes (a) for models SLyl313a and (b) 
SLyl25135a. The time is shown in the upper right corner of each panel in units of ms. fl along x and y axes is computed by 
/x and jy, and hence, it is not determined at the origin. Since the formed hypermassive neutron star is not spheroidal, Q 
along two axes is significantly different. 
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merger sets in. On the other hand, for models with mass 
smaller than a threshold mass iWthr) a hypermassive neu¬ 
tron star is formed at least temporarily. The hypermas¬ 
sive neutron star is stable against gravitational collapse 
for a while after its formation, but it will collapse to a 
black hole eventually due to radiation reaction of grav¬ 
itational waves or due to outward angular momentum 
transfer (see discussion later). 

In the formation of the hypermassive neutron stars, a 
double core structure is first formed, and then, it relaxes 
to a highly nonaxisymmetric ellipsoid (cf. Figs. 0113 and 
ETb)). The contour plots drawn for a high-density region 
with p > 4 X 10^^ g/cm^ show that the axial ratio of the 
bar measured in the equatorial plane is ^ 0.5; the axial 
lengths of the semi major and minor axes are ~ 20 and 

10 km, respectively. FigureEDb) also shows that the axial 
length along the z axis is about 10 km. Namely, a highly 
elliptical rotating ellipsoid is formed. This outcome is 
signihcantly different from the previous ones found with 
the F = 2 equation of state , in which nearly axisym- 
metric spheroidal neutron stars are formed. The reason 
is that the adiabatic index of the realistic equations of 
state adopted in this paper is much larger than 2 that is 
adopted in the previous one. According to a Newtonian 
study 9, uniformly rotating ellipsoid (Jacobi-like el¬ 
lipsoid) exists only for F ^ 2.25. This fact suggests that 
rapidly rotating stars with a large adiabatic index are 
only subject to the ellipsoidal deformation. Note that 
similar results have been already reported in Newtonian 
and post Newtonian simulations [^1^1^ . 

The rotating hypermassive neutron stars also oscillate 
in a quasiradial manner (cf. Fig. (TJ. Such oscillation 
is induced by the approaching velocity at the collision of 
two stars. By the radial motion, shocks are formed at the 
outer region of the hypermassive neutron stars to convert 
the kinetic energy to the thermal energy. The shocks are 
also generated when the spiral arms hit the oscillating 
hypermassive neutron stars. These shocks heat up the 
outer region of the hypermassive neutron stars for many 
times, and as a result, the thermal energy of the envelope 
increases fairly uniformly. The further detail of these 
thermal properties is discussed in Sec. inmai 

Since the degree of the nonaxial symmetry is suffi¬ 
ciently large, the hypermassive neutron star found in 
this paper is a stronger emitter of gravitational waves 
than that found in |l2| . The signihcant radiation de¬ 
creases the angular momentum of the hypermassive neu¬ 
tron stars. The nonaxisymmetric structure also induces 
the angular momentum transfer from the inner region to 
the outer one due to the hydrodynamic interaction. As 
a result of these effects, the rotational angular velocity 

11 = decreases and its profile is modified. In Fig. 
we show the evolution of Q of the hypermassive neutron 
stars along x and y axes at t = 3.897, 6.069, and 8.621 
ms for models SLyl313a and SLyl25135a. At its for¬ 
mation, the hypermassive neutron stars are strongly dif¬ 
ferentially and rapidly rotating. The strong differential 
rotation yields the strong centrifugal force, which plays 


an important role for sustainin g th e large selfgravity of 
the hypermassive neutron stars [23, Ej] . Since the angu¬ 
lar momentum is dissipated by the gravitational radia¬ 
tion and redistributed by the hydrodynamic interaction, 
n decreases significantly in the central region, and hence, 
the steepness of the differential rotation near the center 
decreases with time. This effect eventually induces the 
collapse to a black hole. 

It should be also noted that fl along two axes is signif¬ 
icantly different near the origin. The reason at t = 3.897 
ms is that the formed hypermassive neutron stars have a 
double-core structure (cf. Figs. 2|and 0 and the angu¬ 
lar velocity of the cores are larger than the low density 
region surrounding them. The reason for t ^ 6 ms is 
that the hypermassive neutron star is not a spheroid but 
an ellipsoid of high ellipticity and the angular velocity 
depends strongly on the coordinate i^. 

Figures 0 and 0 show that even after the emission of 
gravitational waves for ~ 10 ms, the hypermassive neu¬ 
tron star is still highly nonaxisymmetric. This indicates 
that gravitational waves will be emitted for much longer 
time scale than 10 ms. Thus, the rotational kinetic en¬ 
ergy and the angular momentum will be subsequently 
dissipated by a large factor, eventually inducing the col¬ 
lapse to a black hole. 

We here estimate the lifetime of the hypermassive neu¬ 
tron stars using Fig. 0 which shows that the value of Oc 
decreases gradually with time. It is reasonable to ex¬ 
pect that the collapse to a black hole sets in when the 
value of ttc becomes smaller than a critical value. Since 
the angular momentum should be sufficiently dissipated 
before the collapse sets in, the threshold value of Oc for 
the onset of the collapse will be approximately equal to 
that of marginally stable spherical stars (i.e., the dot¬ 
ted horizontal line in Fig. 0 . One should keep in mind 
that the threshold value depends on the slicing condition, 
and thus, this method can work only when the same slic¬ 
ing is used for computation of the spherical star and for 
simulation. In this paper, the (approximate) maximal 
slicing is adopted both in the simulation and in compu¬ 
tation of spherical equilibria so that this method can be 
used. The results for models SLyl35135b and FPS1212b 
indeed illustrate that the prediction by this method is 
appropriate. 

For models SLyl313a, SLyl25135a and SLyl212b, the 
decrease rate of the value of Oc estimated from the data 
for 5 ms ^ t ^ 10 ms is ~ 0.005 ms“^. Extrapolating 
this result suggests that the hypermassive neutron stars 
will collapse to a black hole at t ^ 30 ms for models 
SLyl313a and SLyl25135a and at t ~ 50 ms for model 
SLyl212b. These time scales are much shorter than the 
dissipation time scale by viscosities or the redistribution 
time scale of the angular momentum by the magnetic 
field 1^. Therefore, the gravitational radiation or the 
outward angular momentum transfer by the hydrody¬ 
namic interaction plays the most important role in the 
evolution of the hypermassive neutron stars. 

In the prompt formation of a black hole, most of the 
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fluid elements are swallowed into the black hole in 1 ms 
after the merger sets in. Thus, the flnal outcome is a 
system of a rotating black hole and a surrounding disk 
of very small mass (cf. Fig. IHIa)). In Fig. El we plot the 
evolution of the total baryon rest-mass located outside a 
radius r, M*(r), for models SLyl414a and SLyl35145a. 
M*(r) is defined by 


M, 


(r)= f 

J r<r'<rj 


p^Sx', 


(39) 


where rmax is introduced to exclude the contribution from 
the small-density atmosphere. In the present work we 
choose as r^ax = 7.5Mo « 31 km within which the inte¬ 
grated mass of the atmosphere is negligible (< 0 . 01 % of 
the total mass). The results are plotted for t/Mq = 1.5, 
3, and 4.5. Note that the apparent horizon is located 
for r « O.SMq at t « 3.0 ms for models SLyl414a and 
SLyl35145a, and inside the horizon about 99% and 98% 
of the initial mass are enclosed for these cases, respec¬ 
tively. Figure E) indicates that the fluid elements still 
continue to fall into the black hole at the end of the sim¬ 
ulation. This suggests that the flnal disk mass will be 
smaller than 1 % of the total baryon rest-mass. 

In ^3, we found that even the small mass difference 
with Qm 0.9 increases the fraction of disk mass around 
the black hole significantly. However, in the present equa¬ 
tions of state, Qm ~ 0.9 is not small enough to signif¬ 
icantly increase the disk mass. This results from the 
difference in the equations of state. The detailed reason 
is discussed in Sec. iTTmii 

The area of the apparent horizons A is determined in 
the black hole formation cases. We And that 


A 

WttMq 


0.85, 


(40) 


for models SLyl414a and SLyl35145a. Since most of the 
fluid elements are swallowed into the black hole and also 
the energy carried out by gravitational radiation is at 
most ^ O.OlMo (cf. Fig. IH^ . the mass of the formed 
black hole is approximately ^ 0.99Mo. Assuming that 
the area of the apparent horizon is equal to that of the 
event horizon, the nondimensional spin parameter of the 
black hole defined by g = Jeh/AI^^, where Jbh and 
Mbh are the angular momentum and the mass of the 
black hole, are computed from 


A 

IbTrMgjj 


l + (l-g2)i/2 


(41) 


Equation (ITO implies that for A/167rM|jj ~ 0.85, g ~ 

0.7. 

On the other hand, we can estimate the value of g 
in the following manner. As shown in Sec. IIII 01 the 
angular momentum is dissipated by ^ 10-15% by grav¬ 
itational radiation, while the ADM mass decreases by 
~ 1%. As listed in Table II, the initial value of g is 
~ 0.9. Therefore, the value of g in the flnal stage should 


be ~ 0.75-0.8. The values of g derived by two indepen¬ 
dent methods agree with each other within ^ 10 % error. 
This indicates that the location and the area of the black 
holes are determined within ^ 10 % error. 

For g = 0.7-0 .8 and Mbh ~ 2.8Mq, the frequency 
of the quasinormal mode for the black hole oscillation is 
about 6 . 5 - 7 ( 2 . 8 M 0 /Mbh) kHz ^5. This value is rather 
high and far out of the best sensitive frequency range of 
the laser interferometric gravitational wave detectors . 
Thus, in the following, we do not touch on gravitational 
waveforms in the prompt black hole formation. 


2. Threshold mass for black hole formation 

The threshold value of the total ADM mass above 
which a black hole is promptly formed is Mthr ~ 2.7Mq 
for the SLy equation of state and Mthr ~ 2.5Mq for 
the FPS one with Fth = 2. For the SLy case, we And 
that the value does not depend on the mass ratio for 
0.9 ^ Qm < 1- The maximum allowed mass for the sta¬ 
ble and spherical neutron stars is 2 .O 4 M 0 and 1 . 8 OAf 0 for 
the SLy and FPS equations of state, respectively. This 
implies that if the total mass is by ~ 30-40% larger than 
the maximum allowed mass for stable and spherical stars, 
a black hole is promptly formed. In a previous study with 
the F = 2 equation of state m , we found that threshold 
mass is by about 70% larger than the maximum allowed 
mass for stable and spherical neutron stars. Thus, com¬ 
paring the threshold value of the total ADM mass, we 
can say that a black hole is more subject to be formed 
promptly with the realistic equations of state. 

In the maximum mass of differentially rotating 

stars in axisymmetric equilibrium (hereafter Mmax:dif) is 
studied for various equations of state. The authors com¬ 
pare Mniax:dif with the maximum mass of spherical stars 
(hereafter M„iax:sph) for given equations of state. They 
And that the ratio M„iax:dif/Mniax:sph for FPS and APR 
equations of state (APR is similar to SLy equation of 
state) is much smaller than that for P = 2 equation of 
state. Their study is carried out for axisymmetric rotat¬ 
ing stars in equilibrium and with a particular rotational 
law, and hence, their results cannot be simply compared 
with our results obtained for dynamical and nonaxisym- 
metric spacetime. However, their results suggest that 
the merged object may be more susceptible to collapse 
to a black hole with the realistic equations of state. This 
tendency agrees with our conclusion. 

The compactness in each neutron star of no rotation 
in isolation is deflned by C = GMsph/RsphC^ where Mgph 
and Rsph denote the ADM mass and the circumferential 
radius of the spherical stars. For the SLy equation of 
state, C « 0.151, 0.165, 0.172, and 0.178 for M^ph = 

1.2, 1.3, 1.35, and 1 . 4 M 0 , respectively. For FPS one, 
C w 0.162, 0.169, 0.177, and 0.192 for M^ph = 1.2, 1.25, 

1.3, and 1 . 4 M 0 , respectively. This indicates that a black 
hole is promptly formed for C ^ 0.17 after merger of two 
(nearly) identical neutron stars. In the F = 2 equation of 
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state, the threshold value of C is ~ 0.15-0.16 0. Thus, 
comparing the threshold value of the compactness, we 
can say that a black hole is less subject to be formed 
with the realistic equations of state. 

The reason that the threshold mass for the prompt 
black hole formation is smaller with the realistic equa¬ 
tions of state may be mentioned in the following manner: 
In the realistic equations of state, the compactness of 
each neutron star is larger than that with the T = 2 
equation of state for a given mass. Accordingly, for a 
given total mass, the binary system at the onset of the 
merger is more compact. This implies that the angu¬ 
lar momentum is dissipated more before the merger sets 
in with the realistic equations of state. In the case of 
the hypermassive neutron star formation, the centrifugal 
force plays the most important role for sustaining the 
large selfgravity. Thus, the large dissipation of the angu¬ 
lar momentum before the merger helps the prompt black 
hole formation. Therefore, a black hole is more subject 
to be formed in the realistic equations of state. 


3. Thermal properties 


In Fig. Iinia)-(c), we show prohles of e and Eth as 
well as that of p along x and y axes at t = 2.991 ms for 
model SLyl414a, at t = 8.621 ms for model SLyl313a, 
and at t = 8.621 ms for model SLyl25135a, respectively. 
The density contour curves at the corresponding time 
steps are displayed in the last panel of Figs. |3H3 Figure 
cnid) shows the evolution of the total internal energy and 
thermal energy defined by 


t/ = y* p^ed^x, 

(42) 

t/th = y P^^tYid 

(43) 


for models SLyl313a, SLyl25135a, and SLyl212b. Note 
that in the absence of shock heating, e should be equal 
to Ecoid- Thus, Eth denotes the specific thermal energy 
generated by the shock heating. 

First, we focus on the thermal property for models 
SLyl313a and SLyl25135a which are the representative 
models of hypermassive neutron star formation. In these 
cases, the heating is not very important in the central re¬ 
gion. This is reasonable because the shocks generated 
at the collision of two stars are not very strong, and 
thus, the central part of the hypermassive neutron stars 
is formed without experiencing the strong shock heating. 
On the other hand, the shock heating plays an impor¬ 
tant role in the outer region of the hypermassive neutron 
stars and in the surrounding accretion disk since the spi¬ 
ral arms hit the hypermassive neutron stars for many 
times. 

The typical value of Eth is 0.0Ic^-0.02c^. Here, we re¬ 
cover c for making the unit clear. In the following, we 



t (msec) 


FIG. 9: Evolution of M*(r)/M* for r = 1.5, 3, and 4.5Mo for 
models SLyl414a (curves except for the dotted curve) and 
SLyl35145a (dotted curves). Mo and M* denote the ADM 
mass and total baryon rest-mass of binary neutron stars at 
t = Dehnition for Mt,{r) is shown in Eq. 

assume that the components of the hypermassive neu¬ 
tron stars and surrounding disks are neutrons. Then, 
the value of Eth = O.Olc^ implies that the thermal energy 
per nucleon is 

9.37 1 „ ) MeV/nucleon. (44) 

yO.OlG'y 

Since the typical value of Eth is 0.01c^-0.02c^, the typical 
thermal energy is 10-20 MeV. This value agrees approx¬ 
imately with that computed in [4^ 1^ . 

Figure irnT dl shows that the total internal energy and 
thermal energy are relaxed to be 

U ~ 0.1M*c^ ~ 5 X 10*^3 erg, (45) 

Uth ~ 0.025M^c^ ~ 1 X 10^3 erg, (46) 

for both models SLyl313a and SLyl25135a. Thus, the 
thermal energy increases up to ~ 25% of the total inter¬ 
nal energy. We note that these values are approximately 
identical between models SLyl313a and SLyl25135a. 
This implies that the mass ratio of Qm ^ 0.9 does not 
significantly modify the thermal properties of the hyper¬ 
massive neutron stars in the realistic equations of state. 

The region of En, ~ O.Olc^ will be cooled via the emis¬ 
sion of neutrinos |4^E3- According to the emis¬ 

sion rate in the hypermassive neutron star with the aver¬ 
aged value of Eth 10-30 MeV is 10®^-10®^ erg/s. Thus, 
if all the amounts of the thermal energy are assumed 
to be dissipated by the neutrino cooling, the time scale 
for the emission of the neutrinos will be 1-10 s. This is 
much longer than the lifetime of the hypermassive neu¬ 
tron stars ^ 100 ms. Therefore, the cooling does not play 
an important role in their evolution. 
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FIG. 10: Profiles of e (solid curves) and eth = e — Scold (dashed curves) as well as that of p along x and y axes (a) at t = 2.991 
ms for model SLy 1414a, (b) at f = 8.621 ms for model SLyl313a, and (c) at t = 8.621 ms for model SLyl25135a. Note that 
the region of r ^ 2 km for panel (a) is inside the apparent horizon (see Figs. Elandl^a)). (d) Evolution of U and f/th in unit of 
the rest-mass energy M*c^ for models SLyl313a (solid curves), SLyl25135a (dotted curves), and SLyl212b (dashed curves). 


Since the lifetime of the hypermassive neutron stars 
^ 100 ms is nearly equal to the time duration of the 
short gamma-ray bursts [^, it is interesting to ask if 
they could generate the typical energy of the bursts. In 
a model for central engines of the gamma-ray bursts, a 
fireball of the electron-positron pair and photon is pro¬ 
duced by the pair annihilation of the neutrino and an¬ 
tineutrino [i^ l50| . In ^3, Janka and Ruffert estimate 
the efficiency of the annihilation as severalxlO“^ for the 
neutrino luminosity ^ 10^^ erg/s, the mean energy of 
neutrino ^ 10 MeV, and the radius of the hypermassive 
neutron star ~ 10 km (see Eq. (1) of ^3)- This suggests 
that the energy generation rate of the electron-positron 
pair is ~ 10®° erg/s. Since the lifetime of the hypermas¬ 
sive neutron stars is ^ 100 ms, the energy available for 
the fireball will be at most ~ 10^° erg. This value is not 
large enough to explain typical cosmological gamma-ray 
bursts. Furthermore, as Janka and Ruffert found ^3, the 


pair annihilation of the neutrino and antineutrino is most 
efficient in a region near the hypermassive neutron star, 
for which the baryon density is large enough (cf. Fig. 
13b)) to convert the energy of the fireball to the kinetic 
energy of the baryon. Therefore, it is not very likely that 
the hypermassive neutron stars are the central engines of 
the typical short gamma-ray bursts. 

Now, we focus on model SLyl414a in which a black 
hole is promptly formed after the merger. Comparing 
Fig. llOl a) with the last panel of Fig. O the region of high 
thermal energy is located along the spiral arms of the 
accretion disk surrounding the central black hole. (Note 
that the region oir ^ 2 km is inside the apparent horizon, 
and hence, we do not consider such region.) The part of 
the matter in the spiral arms with small orbital radius 
r ^ 5 km is likely to be inside the radius of an innermost 
stable circular orbit around the black hole, and hence, be 
swallowed into the black hole. Otherwise, the matter in 
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the spiral arms will form an accretion disk surrounding 
the black hole. Thus, eventually a hot accretion disk will 
be formed. However, the region of high thermal energy 
for r ^ 10 km is of low density with p ^ 10^^ g/cm^, 
and the total mass of the disk will be IO^^Mq ^ M* ^ 
O.OlMo (see Fig. ©. The total thermal energy of the 
accretion disk is estimated as 

t/th ^ 

where Mdisk and eth denote the mass of the accretion disk 
and the averaged value of the specific thermal energy. 
Hence, even if all the amounts of the thermal energy are 
dissipated by the emission of neutrinos, the total energy 
of the radiated neutrinos will be at most several x 10®° 
erg. According to , the efficiency of the annihilation 
of the neutrino and antineutrino is several x 10“® for the 
neutrino luminosity ^ 10®° erg/s, the mean energy of 
neutrino ~ 10 MeV, and the disk radius ~ 10 km. This 
indicates that the energy of the fireball is at most ^ 10^® 
erg. Although the density of the baryon at the region 
that the pair annihilation is likely to happen is small 
enough to avoid the baryon loading problem, this energy 
is too small to explain cosmological gamma-ray bursts. 


4- Effects of mass difference 

Comparing the evolution of the contour curves, the 
maximum density, and the central value of the lapse func¬ 
tion for models SLyl313a and SLyl25135a (see Figs. 
and0a)). it is found that the mass difference plays a mi¬ 
nor role in the formation of a hypermassive neutron star 
as far as Qm is in the range between 0.9 and 1. Figures 
0 a) and 0 also illustrate that the evolution of the system 
to a black hole is very similar for models SLyl414a and 
SLyl35145a. In in which simulations were performed 
using the F = 2 equation of state, we found that the mass 
difference with Qm ~ 0.9 significantly induces an asym¬ 
metry in the merger which contributes to formation of 
large spiral arms and the outward angular momentum 
transfer, which are not very outstanding in the present 
results. The reason seems to be as follows. In the previ¬ 
ous equation of state, the mass difference with Qm 0.9 
results in a relatively large (^ 15%) difference of the 
compactness between two stars. On the other hand, the 
difference in the compactness between two stars with the 
present equations of state is ~ 10% for Qm ~ 0.9. This 
is due to the fact that the stellar radius depends weakly 
on the mass in the range 0.8M© ^ M ^ I.5M© (see Fig. 
0. As a result of the smaller difference in the compact¬ 
ness, the tidal effect from the more massive star to the 
companion becomes smaller, and therefore, the asymme¬ 
try is suppressed. To yield a system of a black hole and 
a massive disk, smaller mass ratio with Qm < 0.9 will be 
necessary in the realistic equations of state. 


Another possible reason is that neutron stars in the 
realistic equations of state are more compact than those 
in the F = 2 equation of state. Due to this fact, at the 
merger, the system is more compact, and hence, even in 
the formation of the asymmetric spiral arms, they cannot 
spread outward extensively but wind around the formed 
neutron star quickly. Consequently, the mass of the disk 
around the central object is suppressed to be small and 
also the asymmetric density configuration does not be¬ 
come very outstanding. 


5. Dependence of dynamics on the grid size and Fth 
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FIG. 11: The same as Fig. |7|but for models SLyl313a (solid 
curves), SLyl313b (dashed curves), SLyl313c (dotted dashed 
curves), and SLyl313d (long dashed curves). Note that the 
simulation for SLyl313b was stopped at t ~ 7.5 ms to save 
computational time. 


For model SLyl313, we performed simulations chang¬ 
ing the value of Fth with the grid size (377, 377, 189). 
In Fig. CH evolution of ttc and pmax is shown for models 
SLyl313b-SLyl313d. Note that the grid size and grid 
spacing are identical for these models. The results for 
model SLyl313a are shown together for comparison with 
those for model SLyl313b for which the parameters are 
identical but for the grid size. 

By comparing the results for models SLyl313a and 
SLyl313b, the magnitude of the error associated with 
the small size of L is investigated. Figure ITTI shows that 
two results are approximately identical besides a system¬ 
atic phase shift of the oscillation. This shift is caused 
by the inappropriate computation of the radiation reac¬ 
tion in the late inspiral stage for t ^ 2 ms: For model 
SLyl313b, L is smaller, and hence, the radiation reaction 
in the inspiral stage is significantly overestimated to spu¬ 
riously accelerate the merger resulting in the phase shift. 
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However, besides the phase shift, the results are approx¬ 
imately identical. In particular, the results agree well in 
the merging phase. This indicates that even with the 
smaller grid size (377, 377, 189), the formation and evo¬ 
lution of the hypermassive neutron star can be followed 
within a small error. 

Comparison of the results among models SLyl313b- 
SLyl313d tells that for the smaller value of Tth, the max¬ 
imum density (central lapse) of a hypermassive neutron 
star formed during the merger is larger (smaller). This 
is due to the fact that the strength of the shock formed 
at the collision of two stars, which provides the thermal 
energy in the outer region of the formed hypermassive 
neutron stars to expand, is proportional to the value of 
Tth - 1. Figure CH also indicates that the evolution of 
the system does not depend strongly on the value of Tth 
for Tth ^ 1.65. However, for Fth = 1-3, the formed hy¬ 
permassive neutron star is very compact at its birth, and 
hence, collapses to a black hole in a short time scale (at 
t 8 ms) after the angular momentum is dissipated by 
gravitational radiation. This time scale for black hole for¬ 
mation is much shorter than that for models SLyl313b 
and SLyl313d for which it would be ~ 30 ms. This indi¬ 
cates that for small values of Tth — 1 ^ 0.3, the collapse 
to a black hole is significantly enhanced. 

In reality, in the presence of cooling processes such 
as neutrino cooling, the adiabatic index Fth decreases 
effectively. Such cooling mechanism may accelerate the 
formation of a black hole. However, the emission time 
scale of the neutrino is ~ 1-10 s as mentioned in Sec. 
irrmiii Thus, the effect does not seem to be very strong. 


C. Gravitational waveforms 


1. Waveforms in the formation of hypermassive neutron 
stars 


In Figs. cmsi we present gravitational waveforms 
in the formation of hypermassive neutron stars for sev¬ 
eral models. Figure CHa) shows R+ and i?x for model 
SLyl313a. For comparison, gravitational waves com¬ 
puted in terms of a quadrupole formula (H-|_ and Hx) 
defined in Sec. 11 B are shown together by the dashed 
curves. The amplitude of gravitational waves, h, ob¬ 
served at a distance of r along the optimistic direction 
{9 = 0) is written as 


h 


10-22 


i?+,x \/100Mpc\ 
0.31 km y \ r ) 


(48) 


Thus, the maximum amplitude observed along the most 
optimistic direction is ~ 2 x 10“^^ at a distance of 100 
Mpc. 

In the real data analysis of gravitational waves, a 
matched filtering technique 0 is employed. In this 
method, the signal of the identical frequency can be accu¬ 
mulated using appropriate templates, and as a result, the 


effective amplitude increases by a factor of where 
N denotes the number of the cycle of gravitational waves 
for a given frequency. We determine such effective am¬ 
plitude in Sec. IHIG3l fcf. Eq. (I^ i. 

The waveforms shown in Fig. da) are typical ones 
in the formation of a hypermassive neutron star. In the 
early phase (tret ^ 2 ms), gravitational waves associated 
with the inspiral motion are emitted, while for tret ^ 2 
ms, those by the rotating and oscillating hypermassive 
neutron star are emitted. In the following, we focus only 
on the waveforms for tret ^ 2 ms. 

Gravitational waves from the hypermassive neutron 
stars are characterized by quasiperiodic waves for which 
the amplitude and the frequency decrease slowly. The 
amplitude decreases with the ellipticity, which is de¬ 
creased by the effects that the angular momentum de¬ 
creases due to the radiation reaction and is transferred 
from the inner region to the outer one by the hydrody¬ 
namic interaction associated with the nonaxisymmetric 
structure. However, the time scale for the decrease ap¬ 
pears to be much longer than 10 ms as illustrated in Figs. 
dCSl The oscillation frequency varies even more slowly. 
The reason seems to be that the following two effects 
approximately cancel each other: (i) with the decrease 
of the angular momentum of the hypermassive neutron 
stars due to the radiation reaction as well as the angu¬ 
lar momentum transfer by the hydrodynamic interaction 
with outer envelope, the characteristic frequency of the 
figure rotation decreases, while (ii) with the decrease of 
the angular momentum, the centrifugal force is weakened 
to reduce the characteristic radius for a spin-up. (We 
note that the radiation reaction alone may increase the 
frequency of the hgure rotation |5l| . In the hypermassive 
neutron stars formed after the merger, the angular mo¬ 
mentum transfer due to the hydrodynamic interaction is 
likely to play an important role for the decrease of the 
frequency.) 

In gravitational waveforms computed in terms of the 
quadrupole formula (the dashed curves in Fig. I12II . the 
amplitude is systematically underestimated by a factor 
of 30-40%. This value is nearly equal to the magni¬ 
tude of the compactness of the hypermassive neutron 
star, GM/Rc^, where M and R denote the character¬ 
istic mass and radius. Since the quadrupole formula is 
derived ignoring the terms of order GM/Rc^, this mag¬ 
nitude for the error is quite reasonable. In simulations 
with Newtonian, post Newtonian, and approximately rel¬ 
ativistic frameworks, gravitational waves are computed 
in the quadrupole formula (e.g., miillll). The re¬ 
sults here indicate that the amplitudes for quasiperiodic 
gravitational waves from hypermassive neutron stars pre¬ 
sented in those simulations are signihcantly underesti¬ 
mated Although the error in the amplitude is large. 


^ Besides systematic underestimation of the wave amplitude, 
rather quick damping of quasiperiodic gravitational waves is seen 
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FIG. 12: (a) Gravitational waves for model SLyl313a. R+ and i?x (solid curves) and A+ and Ax (dashed curves) as functions 
of the retarded time are shown, (b) R+ and i?x as functions of the retarded time for model SLyl25135a (solid curves). For 
comparison, those for SLyl313a are shown by the dashed curves. 


the wave phase is computed accurately except for a slight 
systematic phase shift. From the point of view of the 
data analysis, the wave phase is the most important in¬ 
formation on gravitational waves. This suggests that a 
quadrupole formula may be a useful tool for comput¬ 
ing approximate gravitational waves. We note that these 
features have been already found for oscillating neutron 
stars |88) and for nonaxisymmetr ic p rotoneutron stars 
formed after stellar core collapse |^. Here, we recon¬ 
firm that the same feature holds for the merger of binary 
neutron stars. 

In Fig. db), we display gravitational waveforms for 
model SLyl25135a. For comparison, those for SLyl313a 
are shown together (dashed curves). It is found that two 
waveforms coincide each other very well. As mentioned in 
Sec. Ill I l3 4l the mass difference with Qm 0.9 does not 
induce any outstanding change in the merger dynamics. 
This fact is also reflected in the gravitational waveforms. 

In Fig. 1131 we compare gravitational waveforms for 
models SLyl313b (solid curves) and SLyl313a (dashed 
curves). For SLyl313b, the simulation was performed 
with a smaller grid size and gravitational waves were ex¬ 
tracted in a near zone with T/Aq « 0.25 and T/Amerger ~ 
0.83 (cf. for model SLyl313a, L/Aq « 0.41 and 
^/Amerger ~ 1.39). This implies that the waveforms for 
model SLyl313b are less accurately computed than those 
for SLyl313a. Indeed, the wave amplitude for tret ^ 2 ms 
is badly overestimated. However, the waveforms from the 


in the results of these references. This quick damping also seems 
to be due to a systematic error, which is likely to result from a 
relatively dissipative numerical method (SPH method) used in 
these references. 
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FIG. 13: R+ and i?x for models SLyl313a (dashed curves) 
and SLyl313b (solid curves) as functions of the retarded 
time. It is found that gravitational waveforms in the merger 
stage agree well. Note that the simulation for SLyl313b was 
stopped at t ~ 7.5 ms to save computational time. 


formed hypermassive neutron stars for two models agree 
very well except for a systematic phase shift, which is 
caused by the overestimation for the radiation reaction 
in the early phase (tret ^ 2 ms). Thus, for computa¬ 
tion of gravitational waves from the hypermassive neu¬ 
tron stars, we may choose the small grid size. Making use 
of this fact, we compare gravitational waveforms among 
several models computed with the small grid size in the 
following. 

In Fig. 1141 we compare gravitational waves from the 
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FIG. 14: R+ and Rx for models SLyl35135b (solid curves) 
and SLyl313a (dashed curves) as functions of the retarded 
time. 


hypermassive neutron stars for models SLyl313a (dashed 
curves) and SLyl35135b (solid curves). As shown in Figs. 
ll2l and [T!^ quasiperiodic waves for which the frequency is 
approximately constant are emitted for model SLyl313a. 
On the other hand, the frequency is not constant but 
modulates with time for model SLyl35135b (e.g., see the 
waveforms at fret ~ 3.6, 4.6, 5.6, and 6.6 ms for which 
the wavelength is relatively short). The reason is that the 
formed hypermassive neutron star quasiradially oscillates 
with a large amplitude and the frequency of gravitational 
waves varies with the change of the characteristic radius. 
Due to this, the Fourier spectra for models SLyl313 and 
SLyl35135 are significantly different although the differ¬ 
ence of the total mass is very small (cf. Fig. EJb)). 

In Fig. CHa), we compare gravitational waveforms 
for models SLyl313b (dotted curves), SLyl313c (solid 
curves), and SLyl313d (dashed curves). For these mod¬ 
els, the cold part of the equation of state is identical but 
the value of Fth is different. As mentioned in Sec. immi 
with the smaller values of Fth, the shock heating is less 
efficient, and as a result, the formed hypermassive neu¬ 
tron star becomes more compact. Since the characteris¬ 
tic radius decreases, the amplitude of gravitational waves 
decreases and the frequency increases. This shows that 
the strength of the shock heating affects the amplitude 
and the characteristic frequency of gravitational waves. 

In Fig. mb), we compare gravitational waveforms for 
models SLyl212b (solid curves) and FPS1212b (dashed 
curves). For these models, the equations of state are dif¬ 
ferent, but the total ADM mass is approximately iden¬ 
tical. Since the FPS equation of state is slightly softer 
than the SLy one, the compactness of each neutron star 
is larger by a factor of 5-10% (cf. Fig. and so is for 
the formed hypermassive neutron star. As a result, the 
frequency of gravitational waves for the FPS equation of 
state is slightly (~ 15%) higher (cf. Fig. [TTI dD. On the 


other hand, the amplitude of gravitational waves is not 
very different. This is due to the fact that with increas¬ 
ing the compactness, the radius of the hypermassive neu¬ 
tron star decreases while the angular velocity increases. 
These two effects approximately cancel each other, and as 
a result, dependence of the amplitude is not remarkable 
between two models. 


2. Emission rate of the energy and the angular momentum 

In Fig. Ilbl a'l. the emission rates of the energy and the 
angular momentum by gravitational radiation are shown 
for models SLyl313a (solid curves) and SLyl25135a 
(dashed curves). In the inspiral phase for t^et ^ 2 ms, 
they increase with time since the amplitude and the fre¬ 
quency of the chirp signal increase. After the peak is 
reached, the emission rates quickly decrease by about 
one order of magnitude since the merged object becomes 
a fairly axisymmetric transient object. However, because 
of its large angular momentum, the formed hypermassive 
neutron star soon changes to a highly ellipsoidal object 
which emits gravitational waves significantly. The lumi¬ 
nosity from the ellipsoidal neutron star is as high as the 
first peak at 2.2 ms. This is in contrast with the 

results obtained with the F = 2 equation of state in which 
the magnitude of the second peak is 30-50% as large as 
that of the first peak . This reflects the fact that the 
degree of the ellipticity of the formed hypermassive neu¬ 
tron star is much higher than that found in because 
of the large adiabatic index for the realistic equations of 
state. 

The emission rates of the energy and the angular mo¬ 
mentum via gravitational waves gradually decrease with 
time, since the degree of the nonaxial symmetry de¬ 
creases. However, the decrease rates are not very large 
and the emission rates at fret ~ 10 ms remain to be as 
high as that in the late inspiral phase as dE/dt ~ 7x 10®"* 
erg/s and dJ/dt ~ 7 x 10®° g cm^/s^. The angular mo¬ 
mentum at t ~ 10 ms is J ~ 0.7Jo ~ 4 x 10^° g cm^/s. 
Assuming that the emission rate of the angular momen¬ 
tum does not change and remains ^ 7 x 10®° g cm^/s, 
the emission time scale is evaluated as J/{dJ/dt) ^ 50 
ms. For more accurate estimation, we should compute 
(J — Jmin)/{dJ/dt) where Jmin denotes the minimum al¬ 
lowed angular momentum for sustaining the hypermas¬ 
sive neutron star. Since Jmin is not clear, we set Jmin = 0. 
Thus, the estimated value presented here is an approxi¬ 
mate upper limit for the emission time scale (see discus¬ 
sion below), and hence, the hypermassive neutron star 
will collapse to a black hole within 50 ms. This estimate 
agrees with the value ^ 30 ms obtained in terms of the 
change rate of ac (cf. Sec. nTnu. Therefore, we con¬ 
clude that the lifetime of the hypermassive neutron stars 
and hence the time duration of the emission of quasiperi¬ 
odic gravitational waves are as short as ~ 30-50 ms for 
models SLyl313a and SLyl25135a. 

Figure irrn' bl displays the emission rates of the energy 
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FIG. 15: R+ and Rx (a) for models SLyl313b (dotted curves), SLyl313c (solid curves), and SLyl313d (dashed curves), and 
(b) for models SLyl212b (solid curves) and FPS1212b (dashed curves) as functions of the retarded time. 
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FIG. 16: dE/dt and dj/dt of gravitational radiation (a) for models SLyl313a (solid curves) and SLyl25135a (dashed curves), 
and (b) for models SLyl212b (solid curves) and FPS1212b (dashed curves). 


and the angular momentum for models SLyl212b (solid 
curves) and FPS1212b (dashed curves). For these mod¬ 
els, the value of L is not large enough to accurately com¬ 
pute gravitational waves in the inspiral phase for tret ^ 2 
ms. Thus, we only present the results for the merger 
phase. The emission rates for SLyl212b are slightly 
smaller than those for model SLyl313a. This results 
from the fact that the total mass of the system is smaller. 
On the other hand, the emission rates for FPS1212b is 
slightly larger than that for SLyl212b. The reason is that 
the FPS equation of state is softer than the SLy one, and 
as a result, the formed hypermassive neutron star is more 
compact and the rotational angular velocity is larger. 


The hypermassive neutron star formed for FPS1212b col¬ 
lapses to a black hole at t ~ 10 ms. This is induced by 
the emission of the angular momentum by gravitational 
waves. However, the collapse time is shorter than the 
emission time scale evaluated by J/{dJ/dt). This is rea¬ 
sonable because the hypermassive neutron star formed 
for model FPS1212 is close to the marginally stable con¬ 
figuration, and hence, a small amount of the dissipation 
leads to the collapse. This illustrates that the time scale 
J/idJ/dt) should be regarded as the approximate upper 
limit for the collapse time scale. 
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3. Fourier power spectrum 


To determine the characteristic frequency of gravita¬ 
tional waves, we carried out the Fourier analysis. In 
Fig. El the power spectrum dE/df is presented (a) for 
models SLyl313a and SLyl25135a, (b) for SLyl313a and 
SLyl35135b, (c) for SLyl313a, SLyl313c, SLyl313d, and 
(d) for SLyl212b and FPS1212b, respectively. Since the 
simulations were started with the initial condition of the 
orbital period ~ 2 ms (i.e., frequency of gravitational 
waves is ~ 1 kHz), the spectrum of inspiraling binary 
neutron stars for / < 1 kHz cannot be taken into ac¬ 
count. Thus, only the spectrum for / ^ 1 kHz should 
be paid attention. In the panel (a), we plot the following 
Fourier power spectrum of two point particles in circular 
orbits in the second post Newtonian approximation |52| : 
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Here, /i and M denote the reduced mass and the total 
mass of the binary, and x = We note that 

the third post Newtonian terms does not signihcantly 
modify the spectrum since their magnitude is ~ O.OI of 
the leading-order term. Thus, the dotted curve should 
be regarded as the plausible Fourier power spectrum for 
/ ^ I kHz. 

Figure ini shows that a sharp characteristic peak is 
present at / = 3-4 kHz irrespective of models in which 
hypermassive neutron stars with a long lifetime 10 
ms) are formed (see also Table HI for the list of the char¬ 
acteristic frequency). This is associated with quasiperi- 
odic gravitational waves emitted by the formed hyper¬ 
massive neutron stars. The amplitude of the peak is 
much higher than that in the F = 2 equation of state 
0 The reason is that with the realistic equations of 
state, the ellipticity of the formed hypermassive neutron 
stars is much larger, and as a result, quasiperiodic grav¬ 
itational waves of a higher amplitude are emitted. Also, 
the elliptic structure of the hypermassive neutron stars is 
preserved for a long time duration. These effects amplify 
the peak amplitude in the Fourier power spectrum. 

The energy power spectra for models SLyl313a and 
SLyI25I35a are very similar reflecting the fact that the 
waveforms for these two models are very similar (Fig. 
m). This indicates that the spectral shape depends 
very weakly on the mass ratio Qm as far as it is in the 
range between 0.9 and 1. On the other hand, three peaks 
are present at / « 2.6, 3.6, and 4.5 kHz in the energy 
power spectrum for model SLyl35I35b (Figs. II 7r hH. 
Thus, the spectral shape is quite different from that for 
model SLyl3I3a although the total mass is only slightly 
different between two models. The reason is that the am¬ 
plitude of the quasiradial oscillation of the hypermassive 
neutron star is very large and the characteristic radius 
varies for a wide range for model SLyI35135b, inducing 


the modulation of the wave frequency. Indeed, the dif¬ 
ference of the frequencies for the peaks is approximately 
equal to that of the quasiradial oscillation ~ 1 kHz. As a 
result, the intensity of the power spectrum is dispersively 
distributed to multi peaks in this case, and the amplitude 
for the major peak at / ^ 3.6 kHz is suppressed. The 
similar feature is also found for models SLyI3I3c and 
FPSI2I2b for which the hypermassive neutron stars col¬ 
lapse to a black hole within ~ 10 ms. 

Figures irTT cl and (d) illustrate that the amplitude and 
the frequency for the peak around / ~ 3-4 kHz depend 
on the total mass, the value of Fth, and the equations of 
state as in the case of gravitational waveforms. Figure 
Ei: c) indicates that for the larger total mass (but with 
M < Mthr), the peak frequency becomes higher. Also, 
with the increase of the value of Fth, the peak frequency 
is decreased since the formed hypermassive neutron star 
becomes less compact. As Fig. HHc) shows, the peak 
frequency is larger for the FPS equation of state than for 
the SLy one for the same value of the total mass. This is 
also due to the fact that the hypermassive neutron star 
in the FPS equation of state is more compact. 

The effective amplitude of gravitational waves ob¬ 
served from the most optimistic direction (which is paral¬ 
lel to the axis of the angular momentum) is proportional 
to yjdE / df in the manner 




^ Vl^+P + l^xP/ 

21 f dE/df 


= 1.8 X 10" 


\10®^ erg/Hz 


1/2 


100 Mpc 
r 



where r denotes the distance from the source, and i?+,x 
are the Fourier spectrum of 7?+,x- Equation im¬ 
plies that the effective amplitude of the peak is about 4-5 
times larger than that at 1 kHz. Furthermore, the am¬ 
plitude of the peak in reality should be larger than that 
presented here, since we stopped simulations at t ~ 10 
ms to save the computational time, and hence, the inte¬ 
gration time ^ 10 ms is much shorter than the realistic 
value. Extrapolating the decrease rate of the angular 
momentum, the hypermassive neutron star will dissipate 
sufficient angular momentum by gravitational radiation 
until a black hole is formed. As indicated in Secs. Illl Mil 
and nil (T21 the lifetime would be ~ 30-50 ms for models 
SLyl313a and SLyl25135a and ~ 50-100 ms for model 
SLyl212b. Thus, we may expect that the emission will 
continue for such time scale and the effective amplitude 
of the peak of / ~ 3-4 kHz would be in reality amplihed 
by a factor of ~ 3^/^-10^/^ « 2-3 to be ~ 3-5 x 10“^^ at 
a distance of 100 Mpc. Although the sensitivity of laser 
interferometric gravitational wave detectors for / > 1 
kHz is limited by the shot noise of the laser, this value is 
larger than the planned noise level of the advanced laser 
interferometer « 10“^^ ®(//1 kHz)^/^ i It will be in¬ 
teresting to search for such quasiperiodic signal of high 
frequency if the chirp signal of gravitational waves from 
inspiraling binary neutron stars of distance r ^ 100 Mpc 
are detected in the near future. 
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FIG. 17: Fourier power spectrum of gravitational waves dE/df (a) for models SLyl313a (solid curve) and SLyl25135a (dashed 
curve), (b) for models SLy 1313a (solid curve) and SLyl35135b (dashed curve), (c) for models SLyl313a (dashed curve), 
SLyl313c (solid curve), and SLyl313d (long-dashed curve), and (d) for models SLy 1212b (solid curve) and FPS1212b (dashed 
curve). Since the simulations are started when the frequency of gravitational waves is ~ 1 kHz, the spectrum for / < 1 kHz 
is not correct. The dotted curve in the panel (a) denotes the analytical result of dE/df in the second post Newtonian and 
point-particle approximation. The real spectrum for / ;^ 1 kHz is approximated by the dotted curves. 


Detection of the quasiperiodic gravitational waves will 
demonstrate that a hypermassive neutron star of a life¬ 
time much longer than 10 ms is formed after the merger. 
Since the total mass of the binary should be determined 
by the data analysis for the chirp signal emitted in the in¬ 
spiral phase , the detection of the quasiperiodic gravi¬ 
tational waves will provide the lower bound of the binary 
mass for the prompt formation of a black hole Mthr- As 
found in this paper, the value of Mthr depends sensi¬ 
tively on the equations of state. Furthermore, the values 
of Mthr (Mthr ~ 2.7Mq and ~ 2.5Mq for the SLy and 
FPS equations of state, respectively) are very close to the 
total mass of the binary neutron stars observed so far . 
Therefore, the merge of mass ^ Mthr is likely to happen 
frequently, and thus, the detection of gravitational waves 
from hypermassive neutron stars will lead to constrain¬ 


ing the equations of state for neutron stars. For example, 
if quasiperiodic gravitational waves are detected from a 
hypermassive neutron star formed after the merger of a 
binary neutron star of mass M = 2.6M©, the FPS equa¬ 
tion of state should be rejected. As this example shows, 
the merit of this method is that only one detection will 
significantly constrain the equations of state. The further 
detail about this method is described in [b^. 


4- Calibration of radiation reaction 

Figure d shows evolution of the ADM mass and the 
angular momentum computed in a finite domain by Eqs. 
(1^ and as well as the violation of the Hamiltonian 
constraint H defined in Eq. for models SLyl313a 
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FIG. 18: Evolution of ADM mass and angular momentum in units of their initial values Mq and Jo, and violation of the 
Hamiltonian constraint (a) for model SLyl313a and (b) for model SLyl414a. In the upper two panels, the solid curves denote 
M/Mo and J/Jo computed from Eqs. and (ESI, while the long dashed curves denote 1 — AE{t)/Mo and 1 — AJ{t)/Jo, 
respectively (see Eqs. 12711 and 128t l. 


and SLyl414a. The solid curves in the upper two panels 
denote M and J while the dashed curves are Mq — AE 
and Jo — A J which are computed from the emitted en¬ 
ergy and angular momentum of gravitational waves. The 
ADM mass and angular momentum computed by two 
methods should be identical because of the presence of 
the conservation laws. The figure indicates that the con¬ 
servation holds within ^ 2% error for the ADM mass 
and angular momentum (except for the case that a black 
hole is present). This implies that radiation reaction of 
gravitational waves is taken into account within ~ 2% in 
our numerical simulation. 

The error in the angular momentum conservation is 
generated mainly in the late inspiral phase with t ^ 2 ms 
in which L is smaller than the wavelength of gravitational 
waves and the radiation reaction cannot be evaluated ac¬ 
curately. To improve the accuracy for the conservation in 
this phase, it is required to take a sufficiently large value 
of L that is larger than the wavelength. On the other 
hand, the magnitude of the error does not change much 
after the formation of the hypermassive neutron stars for 
t ^ 2 ms as found in Fig. EDa). This implies that the 
radiation reaction of gravitational waves to the angular 
momentum for the formed hypermassive neutron stars is 
computed within 1% error. 

The bottom panels show that the violation of the 
Hamiltonian constraint is of order 0.01 in the absence of 
black holes. Also noteworthy is that the violation does 
not grow but remain small in the absence of black holes. 
This strongly indicates that simulations will be continued 
for an arbitrarily long duration for spacetimes of no black 
hole. On the other hand, the computation crashes soon 
after the formation of a black hole for model SLyl414a. 


This is mainly due to the fact that the resolution around 
the black hole is too poor. If one is interested in the 
longterm evolution of the formed black hole, it is ob¬ 
viously necessary to improve the resolution around the 
black hole to overcome this problem. In the current sim¬ 
ulation, the radius of the apparent horizon is covered 
only by ~ 5 grid points. The axis ym metric simulations 
for black hole formation (e.g., |2^|^) have experimen¬ 
tally shown that more than 10 grid points for the radius 
of the apparent horizon will be necessary to follow evo¬ 
lution of the formed black hole for SOM ~ 0.4 ms. To 
perform such a better-resolved and longterm simulation 
with L ^ O.SAo, the grid size more than (1500, 1500, 750) 
is required, implying that a powerful supercomputer, in 
which the memory and the computational speed are by 
a factor of ^ 10 as large as those of the present com¬ 
putational resources, is necessary. If the computational 
resources are not improved in the near future, adopting 
the adaptive mesh refinement technique will be inevitable 
for following the evolution of the black hole . 


IV. SUMMARY 

We performed fully general relativistic simulations for 
the merger of binary neutron stars adopting realistic 
equations of state. Since the stiffness is significantly dif¬ 
ferent from that in the F = 2 equation of state adopted in 
the previous works (e.g., |l3)j several new features have 
emerged. The following is the summary of the results 
obtained in this paper: 

1: If the ADM mass of the system is larger than ~ 2.7Mq 
(~ 2.5Mq), a black hole is promptly formed in the SLy 
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(FPS) equation of state. Otherwise, a hypermassive neu¬ 
tron star of ellipsoidal shape is formed. This indicates 
that the threshold mass depends on the equations of state 
and the values are very close to those for observed binary 
neutron stars. 

2: In the black hole formation, most of mass elements 
are swallowed into the horizon, and hence, the disk mass 
around the black hole is much smaller than 1 % of the 
total baryon rest-mass as far as the mass ratio Qm is 
larger than 0.9. Although the disk is hot with the thermal 
energy ^ 10-20 MeV, the total thermal energy which 
is available for the neutrino emission is expected to be 
at most ^ 10^° erg. Since the pair annihilation of the 
neutrino and antineutrino to the electron-positron pair 
would be < 10 “'^ seems to be very difficult to 

generate cosmological gamma-ray bursts in this system. 
3: The nondimensional angular momentum parameter 
(J/M^) of the formed Kerr black hole is in the range 
between 0.7 and 0.8. Then, for the system of mass ^ 
2 . 8 M 0 , the frequency of gravitational waves associated 
with quasinormal mode ringing of Z = m = 2 modes 
would be ~ 6.5-7 kHz, which is too high for gravitational 
waves to be detected by laser interferometric detectors. 
4: The hypermassive neutron stars formed after the 
merger have a large ellipticity with the axial ratio ^0.5. 
They rotate with the period of ^ 0.5-1 ms, and thus, be¬ 
come strong emitters of quasiperiodic gravitational waves 
of a rather high frequency / ^ 3-4 kHz. Although the 
frequency is far out of the best sensitive frequency range 
of the laser interferometric gravitational wave detectors, 
the effective amplitude of gravitational waves is very high 
as several x 10“^^ at a distance of r ~ 100 Mpc. Thus, if 
the merger happens for r < 100 Mpc, such gravitational 
waves may be detectable by advanced laser interferome¬ 
ters. The detection of these quasiperiodic gravitational 
waves will be used for constraining the equations of state 
for nuclear matter. 

5: Because of the larger emission rate of gravitational 
waves, the angular momentum of the hypermassive neu¬ 
tron star is dissipated in a fairly short time scale ^ 100 
ms for the mass M ^ 2A-2.1Mq. Also, due to a high 


degree of nonaxial symmetry, the angular momentum is 
transferred outward by the hydrodynamic interaction. 
As a result of these effects, the hypermassive neutron 
stars collapse to a black hole within 100 ms. This time 
scale is much shorter than the viscous dissipation time 
scale and the transport time scale of the angular momen¬ 
tum by magnetic fields [^. Therefore, the gravitational 
radiation or the outward angular momentum transfer by 
the hydrodynamic interaction plays the most important 
role. 

6 : The thermal energy of the outer region of the hyper¬ 
massive neutron stars is high as ~ 10-20 MeV, and the 
total emission rate of the neutrino energy is estimated 
as ^ 10®^ erg/s. The thermal energy is generated by 
the shocks due to the multiple collisions between the spi¬ 
ral arms and the oscillating hypermassive neutron star. 
Thus, the hypermassive neutron star will be a strong 
emitter of neutrinos. However, the emission time scale is 
~ 1-10 s which is much shorter than the lifetime < 100 
ms. This implies that the neutrino cooling plays a minor 
role in the evolution of the hypermassive neutron star. 

7: The mass difference with the mass ratio Qm ~ 0.9 
does not modify the dynamics of the merger and the 
outcome after the merger significantly from that with 
Qm = 1- This disagrees with the previous result which 
was obtained in the simulations performed with the T = 2 
equation of state ■ The reason is that with the realis¬ 
tic equations of state, the radius of neutron stars is small 
as ~ 11-12 km depending weakly on the mass in contrast 
to that in the T = 2 equations of state. 
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